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Chapter 1 Introduction and Literature Review 


1.1. General Background. 

The proposed study is an investigation into the nonlinear oscillations of liquid 
bridges. The liquid bridge is a finite length fluid column which is held by both 
wetting forces (capillarity) and surface tension between two solid end disks which in 
this study are taken to be coaxial. 

In addition to begin of fundamental interest, the liquid bridge configuration serves 
to model floating zones in crystal growth applications. Moreover, recent studies of 
Marangoni convection (including numerical studies) have been performed in this con- 
figuration. Although recent attention has been focused on the fundamental issues 
of liquid column formation and the possible paths to breakage (Padday, 1992, Mar- 
tinez, 1984), the technical applications also provide impetus for investigation into the 
nonlinear dynamics of the liquid bridge. 
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The focus-of this investigation is on aspects of the nonlinear behavior of the finite 
length liquid column. The emphasis is on bridge dynamics. This is in contrast to the 
quite interesting studies which look at possible equilibrium shapes and their (static) 
stability, as illustrated by the recent work of Langbein (1992) into the static stability 
of liquid bridges held between parallel plates. Of course, the initial interface shape of 
the liquid bridge in the dynamical studies will have to satisfy the capillary equation 

and be statically stable. 

The proposed investigation will concentrate on nonlinear fluid dynamics. Thermal 
and solutal fields, so necessary to actual crystal growth process (Brown, 1988), will 
be absent. However, the floating zone milieu in which crystal growth would occur 
does involve fluid dynamics; and results of the proposed investigation should yield 
insight into nonlinear effects which would impact crystal growth. 

In the microgravity environment provided by space shuttle, residual accelerations 
which could affect the stability limits of the liquid bridge configuration exist. In one 
space experiment, an amphora type liquid bridge configuration resulted instead of the 
c-mode4ype that the experiment was designed to excite (Martinez, 1984). The inter- 
face response to periodic residual accelerations oriented parallel to the longitudinal 
axis of the bridge has been studied in a series of theoretical investigations (Zhang k 
Alexander, 1990, Lyell, 1991, Meseguer k Perals, 1991). Only the work of Zhang and 
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Alexander incorporated any nonlinear effects; however, their formation used a sim- 
plified one-dimensional slice model for the liquid bridge. With regard to the response 
of the liquid column to accelerations in alternate orientations relative to the column, 
the work of Lyell (1993) investigates interface stability in the presence of periodic 
forcing aligned normal to the longitudinal axis. However, the column was taken to 
be infinite in extent. 

However, before any detailed investigations into the full nonlinear dynamics of the 
liquid bridge (with resonance effects and external forcing) can be undertaken, it is 
necessary to investigate the effect of nonlinearity on the oscillation frequencies of the 
system. 

With regard to liquid bridge dynamics, it is only recently that the natural oscilla- 
tion frequencies of the liquid bridge have been calculated in an inviscid Unear analysis 
(Sanz, 1985 (axisymmetric); Sanz & Lopez-Diez, 1989 (non-axisymmetric)). A pre- 
liminary attempt has been performed by Eidel and Bauer (1988) to try to determine 
the effect of the nonlinearity upon the natural oscillation frequencies of the liquid 
bridge. However, they did not include the anchored triple contact Une boundary con- 
dition 4 b their formulation, and so have not modeled the finite liquid bridge correctly, 
(in essence, their results have appUcation to an infinite column.) 

The task of the proposed investigation is to determine the effect of nonUnearity 
upon the natural oscillation frequencies of the liquid bridge; to ascertain whether 
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the liquicLbridge system exhibits softening (or hardening) oscillations for a range of 
bridge slenderness parameters as well as higher order frequencies. Moreover, results 
of such an analysis will yield additional information on nonlinear corrections to the 
interface shape. 

The liquid bridge configuration has both the fluid interface as well as solid bound- 
aries provided by the end disks. It can be viewed as a configuration lying between 
the containerized fluid for which only the top boundary is that of a free surface (in a 
1 gravity environment) and the free liquid drop, which involves no solid boundaries. 

1.2. Literature Review. 

1.2.1. Fluid Physics Literature. 

Early work on the oscillation, dynamics, and breakage of liquid bridges utilized 
the one-dimensional inviscid slice model, which is valid in the limit of slender liquid 
bridges. Such a model assumes that axial velocity in the bridge is independent of the 
radial coordinate. Oscillatory frequencies for the slender limit case were determined in 
a linearized analysis by Meseguer (1983), in a study which investigated the dynamics 
of slender liquid bridges using the one-dimensional slice model. Additional efforts 
which utilized the one-dimensional slice model include that of Rivas and Meseguer 
(1984) and Meseguer and Sanz (1985). A report on liquid bridge breakage aboard 
Spacelab-Dl was given by Meseguer et al (1986). 
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It has been roughly a decade since the oscillation frequencies of the infinite length 
cylindrical liquid column were determined (Bauer, 1982). Such calculations were 
extended to include the oscillation frequencies of viscoelastic infinite length cylindrical 
columns (Bauer, 1986). The determination of the natural oscillation frequencies of 
the liquid bridge without the restriction to the slender limit has been performed 
more recently. A three dimensional linear model was developed by Sanz (19S5); 
axisymmetric modes (Sanz 1985) and non-axisymmetric modes (Sanz and Lopez- 
Diez, 1989) have been investigated and oscillation frequencies determined. 

In addition to fundamental interest, information regarding oscillation models is 
important for applications such as crystal growth via float-zone process. Also, recent 
experiments and numerical studies have investigated thermocapillary flow in the liquid 
bridge configuration. (See, for example, Preisesser et al, 1990, and Velten et al, 1991). 

An attempt to evaluate the nonlinear effects on the frequencies of oscillation of 
a liquid bridge has been made by Eidel and Bauer (1988). However, their analysis 
did not impose the restriction of an anchored triple contact line at end disks, and so 
did not characterize correctly the finite length fluid bridge. (In essence, their analysis 
refers to the infinite length cylindrical fluid column). 

It is this problem which the proposed effort will solve. Nonlinear corrections to 
the oscillation frequencies will be determined over a wide range of slenderness param- 
eters (no restriction to the slender limit). Whether the system exhibits hardening or 
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softening characteristics will be ascertained. The nonlinear correction to the interface 
shape will be found. The proposed methodology will utilize a Lindstedt- Poincare ex- 
pansion of the frequency coupled with domain perturbation techniques. 

1.2.2. Mathematical Formulation Literature Review. 

The investigation of the nonlinear corrections to the frequencies of oscillation as 
well as to the interface shape will proceed via utilization of a Lindstedt-Poincare 
expansion for the frequency, and domain perturbation techniques. As the solution 
is required to be periodic in time, the Lindstedt-Poincare technique is applicable. 
Since the domain in which the solution is to be obtained is nal stationary, domain 
perturbation techniques are quite useful. 

Domain perturbation techniques were developed in the context of non-linear water 
wave theory. Early contributions include the work of Tadjbakhsh and Keller (1960), 
and Verma and Keller (1963) in their investigations of standing finite amplitude 
surface waves (in two and three dimensions, respectively). A formal discussion of the 
methodology awaited the contribution of Joseph (1973). However, it was not until the 
1982 arlicle by Lebovitz that a formal equivalence between expansions produced via 
domain perturbation techniques and those derived from more classical techniques (see 
Wehausen and Laitone, 1960) was exhibited. Recently Gu et al. (1988) used domain 
perturbation techniques in their study of the resonant surface waves in a rectangular 
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container subject to periodic forcing oriented in the vertical direction. 

For a further discussion on the genesis and details of the method, see the afore- 
mentioned references. Application of the domain perturbation techniques results in 
the development of a hierarchy of equations at successive orders in a given expan- 
sion parameter. Each of these sets of equations (including governing equations and 
boundary/interface conditions) is to be solved on a fixed domain. This fixed domain 
is achieved via a transformation of the oscillating interface boundary. 

1.3. Objectives. 

The primary objectives of this work are: 

(1) to determine the nonlinear corrections to the interface shape of a naturally 
oscillating finite length liquid column, and 

(2) to determine the nonlinear corrections to the oscillation frequencies for various 
modes of oscillation. The modes of oscillations themselves may be quickly character- 
ized physically by the number of half-waves present upon the free surface. 

The work will investigate the nonlinear characteristics of free oscillations only. 
This not only a very demanding task, it is also the first task which must be accom- 
plished in any rational plan of investigation into the nonlinear behavior of the liquid 
bridge. 

In order to accomplish the objective of the proposed investigation, several subtasks 
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become critical to the overall effort. First among these is the selection of a methodol- 
ogy which may be applied to the governing equations and allow for the incorporation 
of nonlinear efforts. Moreover, the methodology should be such that known linear 
results may also be recovered. 

It is planned that the approach be theoretical and analytical (as opposed to com- 
putational). A methodology capable of achieving the objectives has been selected. It 
is discussed in Chapter II. 

Application of the methodology (Lindstedt-Poincare expansion in conjunction 
with the domain perturbation technique) results in an hierarchical system of equa- 
tions. The system discussed in Chapter III represents a recovery of known Unear 
results. 

Nonlinear corrections to the interface shape are achieved by solution of the second 
order system given in Chapter IV. Graphical results are presented in Chapter V. 

In order to ascertain the nonUnear correction to the natural frequencies of os- 
cillation and thereby satisfy the final objective of this proposed investigation, it is 
necessary to develop a solvability condition at third order in the hierarchical sets of 
equations. Theoretical results are presented in Chapter VI. PreUminary numerical 
results are given in Chapter VII. 


8 



Chapter 2 Governing Equations and 

Formulations 


2.1. Governing Equations. 

The time-periodic, irrotational and incompressible motion of an inviscid fluid col- 
umn of finite length is considered. The following nondimensionalized quantities are 
used: 



2C — spatial coordinates 

U_ — ■ velocity field 

P — pressure 

t — time 


9 



u — angular frequency 

R — radius of the column 

a — interface tension 

p — density of the fluid comprising column 


and where the tildes indicate corresponding dimensional quantities. 

The volume of the undisturbed column is V = irR 2 L, where L is the length of the 
column. The surface of the column during axisymmetric oscillations is described by 
RF{z, f), where F(z,t) is the dimensionless shape function of the column. 

The nondimensional slenderness of the column is defined as 



Using the nondimensionalizations, the equations governing the inviscid time peri- 
odic motion are listed below, where $(r, z, <) denotes the velocity potential function. 

V 2 $(r,z,f) = 0, (0 < r < F(z,t ), -A < z < A), (2.1) 


which results from using the potential form of the velocity field in the conservation 
of mass equation. 

The conservation of momentum equation is given by 



-VP, (0 <r<F, -A < z < A). 


( 2 . 2 ) 
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The hpundary/interface conditions are 


— = 0, (r = 0, -A < z < A). 


(2.3a) 


dF{z,t) d$ d$ dF n . JU 

-“ft 2 T~z~d~z = °’ (r = F( ‘ M)) - 


(2.36) 


A P = 


F 


1 + 


dF_' 

dz . 


1 + 


m 


d 2 F 

' dz 2 

I , ( r = F(z, t)). 


(2.3c) 


V$(r, z, t + 2 rr) = V$(r, z, t). 


(2.3d) 


5$ 

dz 


= 0 . 


t=±A 


(2.3c) 


/ A F 2 (z,f)d2 = 2A. 

J *— A 


(2.3/) 


F(±A,t) = l. (2.3$) 

Equation (2.1) is the Laplace equation governing the flow; (2.2) is the Euler equa- 
tion; (2.3a) is the condition for a zero radial velocity at the center of the column, 
required for the restriction to axisymmetry; (2.3b) and (2.3c) are the kinematic con- 
dition and the normal force equation, respectively, at the interface; (2.3d) is the 
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condition for: periodicity in time of the velocity field; (2.3e) is the condition for zero 
normal velocity at the end disks; ( 2 . 3 f ) is the conservation of the volume condition 
and (2.3g) is the anchored triple contact line condition. 

2.2. Linstedt-Poincare Expansion with Domain Perturbation Method. 

The unknowns of the equations (2.1) - (2.3) are the shape function F(z,t), the 
velocity potential function $(r, z,t), and the frequency u>. These variables will be 
calculated as the terms in expansion of the amplitude of the motion by the Linstedt- 
Poincare method (see Nayfeh & Mook 1979). The dependence of the velocity potential 
$(r, z,t) on the shape of the mathematical domain which is given by the moving 
boundary r = F(z, t) is very complicated. The domain perturbation technique as 
detailed by Joseph (1973) will be applied. 

To immobilize the boundary shape, introduce the change of variables r = fiF(z, t). 
Let e denote a small positive real number. The expansions of the dependent variables 
can be written in terms of e as follows: 
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\ 

F(z,t;e) 


F^{z,t) 

$(r,z,t:e) 

oc k 

* W (/|,*,0 

P(r,z,t;e) 

Jk = 0 ** 

pW(/i,z,0 

K w{t) 1 


V / 


where 


f f '*»(*, 0 


dt k 


P [k \n,z,t) 


d k ${n,z,t\ 0) 
** 

d fc P(/i,z,<;0) 

dc* 


W w = 

- ** 


The static cylindrical column is recovered as the zeroth order solution of the equation, 


with: 


F<°>(z, 0 = 1, * (0) (M, *, 0 = 0, P {0) (n, z, t) = 0. (2.5) 


Using the chain rule for differentiation, each term z,0 and Pl k \fi,z,t) in the 

expansion for the potential and the pressure can be written as a sum of contributions 

dF 

evaluated on the cylindrical domain (0 < ft < 1), and (—A < z < A). Let = 
F^ l l(z,0‘ Then the first few relationships are: 


$ l °l(/i,z,*;0) = $ (o) (/i,z,0 


13 



<*!*)(//, M;0) = «I»< 1 >(/i,2,0 + F< 1 >(2,<) 


~di 




d$ (,) 

dn 


* [3 W.*;0) = + + 

on on 


, v d 3 &0) , . , „ 5 2 *(°) 

+ 3(F< , >(z,0) 1 ^W+3F< 1 >(M)f (2, (M)- 


d/i 3 


dn 1 


... ,,, 3$(®) 

+ ) <*.0) 35 TTi- + f’ l3 w)- 


d/1 3 


dp 


Similarly, the following results are obtained. 


P®(?,z, P,0) 3 P m (ii,z,t) 

/)D(0) 

pW(^,z,i;0) = P<»(r, ,,t) + F*%,t)-f- 

8 pi 1 ) 

P |2 '((i,i,l;0) = P' 1 >(^,r,i) + 2F I "( J ,() £ 3— 


: + F (3 ^,0^ + (F ( "( 2 ,<)) 3 ^ 

P M (<*,*>‘;0) = P m ^,z,t) + 3F"\z,t)^- + 3F^(z,t)^- 


+ 3 (/•<»(*, i)) 3 ^^ + 3f«(z,l)F (2, (z, 
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+ (F^(zJ)) 3 ^- + F^(z,t) dPi ° ] 


<V 




d k $ ftP 

where <^*'(/i,z,f) = -r-r. P^ k \n, z,t) = — — are always defined in the cylindrical 

at* oe* 


coordinates system (0 < /i < 1, —A < z < A). 


2.3. Hierarchical Systems of Equations. 

In this section, the complete system of equations occurring at each order in c are 
displayed. 


Linear System. 

For 0(c), which represents the linear order, we obtain the following governing 
system of equations: 

The Laplace equation: 


V 2 $ w (n,z,t) = 0, (0 < fi < 1, -A < 2 < A). 


Condition on radial velocity required for axisymmetry: 


5/i 


= 0, (/i = 0, -A < z < A). 


Kinematic condition: 


, 0 )5F< 1 > n , 

-9T + -ar= 0 ' (,, = 1) - 
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Normal force balance: 




Periodicity in time: 


V$ (1) (/i,z,< + 2tt) = 


Zero normal velocity required at end disks: 


d$ (1) 


dz 


= 0 . 


:=±A 


Conservation of volume: 


/* FW(z,t)dz = 0. 

J — A 


Triple contact line condition: 


F (1) (±A,0 = 0. 


System at Nonlinear Order c 2 . 

ForO(c J ), which involves the first nonlinear contributions, the following governing 
system of equations has been obtained: 

The Laplace equation: 

V J $ (2) 0i,2,t) = 0, (0 < n < 1, -A < z < A). 
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Condition on radial velocity required for axisymmetry: 


d<P (2) 

dn 


= 0, (fi = 0, -A < 2 < A). 


Kinematic condition: 


(0] dFW d& 2 > 
-uA 0) — — + — — 
dt dfi 


n d)d F(l] n d$WdFW , 

= 2uA l) — — 2F (l) — — 1- 2— — , (n 


dt 


Normal force balance: 


u> 


(o) 


d$ (2) 

nr 


_ fW _ 


<v 


d 2 F (2) 

nr~ 


dz dz 


= 1 ). 


.2u,<»— - 2w<°> - ( f^y 

dt dpdt \ dfi J \ dz ) 


-2(F<") J +^^j +BJV' 1 ', (,1 = 1). 


Periodicity in time: 


V$ (2) (/i, z, t + 2tt) = V<& (2) (/<, *, t). 


Zero normal velocity required at end disks: 

d$ (2) 


dz 


= 0. 


r=±A 


Conservation of volume: 


/ A [*<*> + (F^y 

«/-A l 


dz = o. 
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Triple contact line condition: 


f (2) (±a,o = o. 


Nonlinear System at Order c 3 . 

For 0(t 3 ), we obtain the following governing system of equations: 
The Laplace equation: 

V 2 $ (3) (/i,z,t) = 0, (0 < p < 1, -A < z < A). 


Condition on radial velocity required for axisymmetry: 


a$( 3 ) 

dfi 


= 0, (p = 0, —A < z < A). 


Kinematic condition: 


dt 


$$( 3 ) 


= 3a,< 2 > 


dt 


+ 3u> (1) 


dF™ 

dt 


- 3 F (,) 


“d^ 2 " 


-3F< 2) 


d 2 $< x > 

dp} 


- 3(F (1) ) 2 


dp 3 


+ 3 


d$ (2 > dF™ 

dz dz 


+6 F (1) 


3F (1) d 2 «I» (1) 

dz dpdz 


dFW d*™ 

+ 3 dz dz ’ 


(**-!). 


18 



Normal force balance: 

,d$< 3) 


u> 


(o>: 


dt 


_ pW - 


d 2 F { 3) 

d? 




dt 


dt 


d^idt 


. 3w (o) f (0^!^ _ 3J°)F(^ - 3 J 0) (F (,) ) 2 ^ — 


-3 


d^idt 

^$( l ) 


dfidt 


dfFdt 


( n d$ (I) W 1} _ 3 a$ (1) d$ (2) 


dn dn ~~ dfi dn 2 
^$(i) 0 2 $O) 


dz dz 


-6Fl l) d * W - 6F (1) F (2) + 6(F (1) ) 3 - 3F (1) 

az a/i3z \ az J 


«S-S¥-*. >-->■ 


Periodicity in time: 


V$ (3) (/i,z,t + 27 t) = V$ (3) (/i,^,0* 


Zero normal velocity required at end disks: 

0$(3) 


dz 


= 0. 


i=± A 


Conservation of volume: 

J S JI/r(3) + f( 1 )F (2) ] dz = 0. 

Triple contact line condition: 


F (3) (±A,0 = 0. 
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Chapter 3 Solutions for the Linear Order 


In this chapter, the linearized problem is solved. Functional results from the linear 
order will appear as nonlinear forcing terms in the higher order systems. The major 
results were first proved by Sanz (1985). In this chapter, Sanz’s results are recovered 
in the notation of the present effort. 

In the previous chapter, the governing equations for the linear order 0(e) were 
listed. They are repeated here: 

The Laplace equation: 

V 2 $ (1) 0<, z, t) s 0, (0 < n < 1, -A < r < A). (3.1) 

Condition on radial velocity required for axisymmetry: 

-5— = 0, (fi = 0, -A < Z < A). (3.2a) 
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— 

Kinematic condition: 


— 

- at + 0 „ -°- '" = ■)' 

(3.2 6) 

— 

Normal force balance: 


— 

* a, - f "’- fli i 

(3-2c) 

— 

Periodicity in time: 


— 

V*% ) r,i + 2 J r) = V$< 1 )( / i,z 1 i). 

(3.M) 

— 

Zero normal velocity at end disks: 


— 

dz ±k ~ °‘ 

(3.2 e) 

— 

Conservation of volume: 


— 

J^F"\z,t)dz = 0. 

(32/) 



Triple contact line condition: 




F<‘)(±A,0 = 0. 

(3.2 S ) 


~ d 

For axisymmetric problems, ^ = 0, and using cylindrical coordinates, Laplace 


equation becomes 



V J $ (1) = ff *— ) * (1) + — $<*> 

fidfi \dfi) dz 2 

(3.3) 

— 
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WWte = T{t)R(p)Z(z) and substitute it into the Laplace 

equation to obtain 

— 

1 d ( 0R\ x „ 

xr -° 

(3.4) 

— 

Z" + XZ =0 

(3.5) 


where A is a constant. 



1 sing the change of variables a = \f\ and £ = a// in ( 

3.4). one obtains the 

— 

modified Basel’s equation: 


— 

2 d 2 R dR , 

" R=0 - 

(3.6) 

— 

Thus for A > 0, the general solution for (3.6) is 


— 

R = AI „(( ) + B *-„((), 


— 

where Iq and Kq are the modified Bessel’s functions of the zeroth order (indicated by 
the subscript) of the first kind. Since the z-axis (r = 0) is part of the domain, B = 0, 


in order to preserve finiteness, and so 


— 

R = AIo(() = AI 0 (an). 

(3.7) 

— 

From (3.5), it is clear that 


— 

Z = E cos(n/X z + A) = E cos(az + A). 

(3.8) 

— 

Thus the solution of is 


— 

= /4/ 0 (a/i)cos(orz + A )T(t). 

(3.9) 

— 
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By the condition of zero normal velocity at the end disks, 

d$ (1) 


dz 


= AI 0 {a t i)(-a)T{t) sin(az + A) 


: = ±A 


nn 


= 0 , 


z±A 


which implies that A = — , for some integer n. 


nn 


Let q„ = /„ = — . It follows 

L A 


® (l ' = Y. A„Io{Ut‘)c°sU(z + A)T((). 


( 3 . 10 ) 


n=0 


Note that Iq(£) = Therefore, 

0$(i) 


dn 


= £ ^nUl(U) Ini* + A)r(t)- 


n=0 


As /j(0) = 0, 




dn 


= 0 , 


M=0 


and so the radial velocity condition is satisfied. Since sin(2/ n A) = sin(rnr) = 0, 
the condition of normal velocity on the end disks is also satisfied. Therefore, the 
spatial form of the solution of has been obtained without knowing the specific 
/-dependence. 

The initial potential is assumed to be zero. Using the periodicity in time condition, 
one can define the time dependence of $U) as gj n t. Hence 

= sin/ -An «>**»(* + A). (311) 

n=0 

By the kinematic condition and normal force balance equations, 


(0) dF™ A , 

-w (0) -^r— + — — = 0, (/i = 1), 


dt 


dn 


(3.12) 
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and 


u> 


(o) 


— ir- - f '" 1 - - B,v m = o (,i = i). 


dt ' dz 2 

Assume that the solution of F (1) is in the form of 


F (1) = Q n {z) cost. 


Substituting (3.14) in the kinematic condition (3.12) yields 


u> (0) Q„(*) + £ A n I' 0 (l n )l n cos l n (z + A) = 0. 


n=0 


Substituting (3.14) in the normal force balance equation (3.13), we have 


<KW + C.M = ^ l01 £ A W /»(/„) cos l„(z + A) - 


n=0 


Solving (3.16) yields 


~ £j(°) 


Q n = a w cos z + 0 {l) sin z - BN (l) + j^A n I 0 (l n ) cos l n (z + A) 

n=0 2 *n 


To study the two expressions of Q n {z) in (3.15) and in (3.17), expand 


cos z = ^ C n cos l n (z + A) 

nsO 


to obtain 


OO 

cos ^ — y~) 

k=i 


2 sin A 


IA(1 -ll) 


COS /j*(z + A) 


sin A 


Similarly, 


^ cos A , / i * \ 

sm 2 = 2w T74 n \ cos *3*-i(* + A). 


t=i A(1 — Qic-i) 


(3.13) 


(3.14) 


(3.15) 


(3.16) 


(3.17) 


(3.18) 

(3.19) 
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Therefore (3.17)- becomes 


m sm A m 2sin A 1 . . . \ 

Qn = -W l) — — jj- cos / 2 jfc(2 + A) 

A A fe = 1 1 “ f 2Jc 


+ £ 1 cosM,(* + A) - £*<■> 

A *Tl 1 “ *2*-l 


00 oj(°) 

+ SZ 1 _ /3 A n I 0 (l n ) cos / n (z 4- A). (3.20) 

n= 0 A ‘n 

For n = 0, (z-independent terms), that /i(0) = 0 and /o(0) = 1 give 

0 = Q (i)£IHA + - BN {1) . 

A 

Physically, BN adjusts the pressure. No adjustment at this linear order is needed 
for physical consistency. Therefore, the value of BN ^ is selected to be zero, and so 

O^aW^+uWAo. (3.21) 

A 

Similarly, for n = 2k > 0, 

An = -2a<V*>^[( W <0 >) 2 /„(/ M ) + Ml - llh)Wn)]~'\ (3.22) 

and for n = 2k — 1 > 0, 

A„_, = -2^'»u,( < »^((u,l 0 )) 1 /„((„.,) + M,(l - /!»_, «('»-.))■'• (3-23) 

This gives the solution of F^(z, t) = cost Q n ( z ) as follows: 

F* l \z, t) = cos 1 1 -a* 1 * 

o° (0) 'l 

+a^ cosz + sinz + £ - A n Io(l n ) cos M 2 + A) > . (3.24) 
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(3.25) 


• . T '** " 

with or 

ros t 00 

Fi'\z, o = ii,E ur 0 (i.) cos /„(* + a). 

W n=l 

Substitute (3.24) into the condition of conservation of volume (3.2f) to obtain 

/ A F^(z,t)dz 
J - A 


= COS 


f m 2 sin A ^ 1 . , . v 

T ~S(i^r smWl+A) 


-A 






sin/ 2 *_i(z + A) 


-A 


00 ( 0 ) 1 

+ X) 1 — » sin + a) 

n=l 1 “ ‘n 


-A 


= 0 

and so (3.24) satisfies the conservation of mass for an incompressible fluid column. 
Substituting (3.24) into the triple contact line condition (3.2g), for z = A, the 

anchored triple contact line requirement yields 

,sin A 


F (1) (A,t) = cost Q n ( A) = cost 


o° (0) 1 

+a (1) cos A + /? (1) sin A + £ TTp AnIt o(U cos /„(2A) ^ = 0; 

n=l * ‘n J 

and for z — —A, 

F* x) (-A,t) = cost Q n (— A) = cost — 
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-2. w<°> 


+ a-^cosA — /? (l) sin A 4- £ ——— A n I 0 (l n )cos{0)) = 0. 

n = l ^ *n J 


(3.27) 


Add (3.26) and (3.27), and use the fact that cos{nir) = ( — l) n to obtain 


+ 0 ei cos a + £ -^-yrAuWn) = 0 


sin A 


2. M<°> 


Substitute (3.22) into (3.28) to obtain 


A — tan A ^ 


k= 1 x ‘2* 


( m <°>) 2 


2 tan A 


*=i 


(1 - Ilk) 


+ Mi - 4) 


= 0 


Similarly, by subtracting (3.27) from (3.26), and by using (3.23), 


\ A tan A + £ 

* =1 (1 -llk-x) 


(w<°>) J 


/. .(0)\2 .1 /i i2 \Ioihk-l) 

r 1 + '“- ,<1 "“- ) wmT) 

Equations (3.29) and (3.30) represent the dispersion relations. 


= 0 


(3.28) 


(3.29) 


(3.30) 


Conclusions for the Solutions of the Order 0(c). For the order O(c), 

OC 

$ (1) (p,z,t) = sint £ A n I 0 {l n n) cos l n (z + A) 

n=0 

and 

{ • \ oo (0) 

_ a (l)!!£_ + Q (l) COS 2 + /3 (1) sin 2 + ^ , 2 A n /o(/ n ) COS l n (z + A) 

A 1 - 'n 

or 

<*"(*.<) = ^E/t n u;((„)cosi„( 2 + A). 


27 



where 


O = a'"^+«< 0 M 0 , 

A 


, sin A 


*»* = ~2a^JV — [(J°') 2 I 0 (l 2k ) + Ml - M(M]-\ 


Au-i = -20 {l) J o) ^{(J o) ) 2 I o (hk-x) 
+ Mi(i -ilk-Mhk-iT 1 -, 


and, setting u > = u/(°) to emphasize the existence of multiple modes, 

( 4°>) 2 


A — tan A ^ 


2 tan A 


‘ sl (l ~l\k) 




= 0, 


-A tan A + £ 


K») a 


Jb=l 


(i - c.) 


(u-ri 1 

■*0t‘2fc— 1 ) 


= 0 . 


Numerical solutions for given A can be obtained from these dispersion (eigenvalue) 
equations. The p = 2 mode can be obtained from (3.30) and the p = 3 mode can be 
obtained from (3.29), where p is the number of half-waves in the interface deformation. 
The modes with an odd (even) number of surface deformations are obtained by using 
equations (3.29) ((3.30), respectively). A root finding technique was used in order 
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to determine the wj, 0) values for each A. The program is listed in the Appendix. 
Graphical results are presented in Figure 3.1 and the numerical results are listed in 
Table 3.1. 

The linear solutions of this order (0(e)) recovers A. Sanz’s solutions (1985). 
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ROOT DETERMINATION BY BISECTION METHOD TO FIND OMEGA 


LAMBDA | OMEGA (P=2 MODE) | OMEGA (P=3 MODE ) 


1.0000 

3.264606138984453400 

7.523140628393082400 

1.2000 

2.290288785914530890 

5.517359210610919050 

1.4000 

1.670226411974154330 

4.203184144471603600 

1.6000 

1.251121188390937440 

3.290224461262170050 

1.8000 

0.953872402882467108 

2.628156961557847950 

2.0000 

0.734523546383059167 

2.131828795651618870 

2.2000 

0.566652921197410286 

1.749636681888166120 

2.4000 

0.433454446539947594 

1.448627121229342450 

2.6000 

0.323269135939490443 

1.206860255054819530 

2.8000 

0.226462506734397689 

1.009208756750524130 

3.0000 

0.129816430336827382 

0.844915285737982580 

3.2000 


0.706097678690390085 

3.4000 


0.586748700766904663 

3.6000 


0.482135179607233483 

3.8000 


0.388141417997999019 

4.0000 


0.300614117118225790 

LAMBDA | 

| OMEGA (P=4 MODE) | 

| OMEGA (P*5 MODE ) | 

1.0000 

12.797858131049062900 

18.777629760152812800 

1.2000 

9.527490411879419700 

14.084428617905016900 

1.4000 

7.383664680146235560 

11.005851738380386000 

1.6000 

5.889552334017784220 

8.858385874800074330 

1.8000 

4.799506245931250750 

7 . 290301709168400810 

2.0000 

3. 976672194738063660 

6.104153871476119210 

2.2000 

3.338290221303994660 

5.181477567644345860 

2.4000 

2.831998994119034000 

4.447332002919017310 

2.6000 

2.423088801453758820 

3.852210648810926230 

2.8000 

2.087696090144671060 

3.362202349416284710 

3.0000 

1.808906068875178570 

2.953358820877600270 

3.2000 

1.574420289854091330 

2.608325526860350460 

3.4000 

1.375094424256466660 

2.314158395188188780 

3.6000 

1.204001042690194370 

2.061181711922818180 

3.8000 

1.055797902929259900 

1.841870305525759920 

4.0000 

0.926300795696662571 

1.650362380063753460 
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LAMBDA OMEGA (P=6 MODE) | OMEGA (P=7 MODE ) 


1.0000 

1.2000 

1.4000 
1.6000 
1.8000 
2.0000 
2.2000 

2.4000 
2.6000 
2.8000 
3.0000 
3.2000 

3.4000 
3.6000 
3.8000 
4.0000 


25.467738030123200600 
19.168957516886681700 
15.039590677929982800 
12 . 161136645430588000 
10.059471455407461300 
8.470268297490264330 
7.233823344526586930 
6.249413915126632580 
5.450582931543062730 
4.791893868648658740 
4.241331957985241540 
3.775746660722400170 
3.378012864456509720 
3.035212736699548720 
2.737425280715521710 
2.476916703999844090 


32.754142031693433500 
24 .712642351512379000 
19.439369813704285400 
15.762671957877159900 
13 . 079191696480551200 
11.049919726683331300 
9.471294339905082180 
8.214513832394541030 
7.194550792016014110 
6.353266606388104480 
5.649717318975978400 
5.054338388853206740 
4.545165435275326570 
4.105857244407967070 
3.723756342366369540 
3 . 389043470572185690 


LAMBDA | 

OMEGA (P«8 MODE) 

1.0000 

40.629465277129657600 

1.2000 

30.697886249731954700 

1.4000 

24.185997431375284300 

1.6000 

19.646925116789443400 

1.8000 

16.333418982145715900 

2.0000 

13.828645817765947300 

2.2000 

11.880501548046547500 

2.4000 

10.329839411058587800 

2.6000 

’ 9.071571819568198690 

2.8000 

8.033823073607900160 

3.0000 

7.165957682358648610 

3.2000 

6.431391677787967830 

3.4000 

5.803106046519861620 

3.6000 

5.260761645784036710 

3.8000 

4.788783435320783880 

4.0000 

4.375064459219278220 
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Chapter 4 Solutions for the Order of 0(e 2 ) 


In this chapter, we shall discuss the solution to the system at the order 0(e 2 ). 
Recall that in a previous chapter, the governing equations for the order 0(( 2 ) were 
listed, and they are now repeated here: 

The Laplace equation: 


V 2 $ (a) (/*»z,t) = 0, (0 < n < 1, -A < z < A). 


(4.1) 


Condition on radial velocity required for axisymmetry: 


d\i 


= 0, (/i = 0, -A < z < A). 


(4.2a) 


Kinematic condition: 


Normal force balance: 


u> 


( 0 ) 


~dT 


_ f (2) - 


d 2 F {2) 

~d?~ 
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= -20^ 


d$ (1) 

~aT' 


2 u {0) F {1) 


d 2 $ {l) 

dfidt 



-2(F(1,)J + (^lr) + 5 - v<2) > (/* = o- 

Periodicity in time: 


v$ (2) (/i, 2 , t + 2tt) = V<J> (2) (//, 2 , 0- 


(4.2c) 


(4-2 d) 


Zero normal velocity at end disks: 


dz 


= 0 . 


(4.2e) 


2=±A 


Conservation of mass condition: 


/« [ F<11 + ( f<1, )1 Jz = °- 


(4.2/) 


Triple contact line condition: 


f ( 2 ) (±a,<) = o. 




In general, a solution of $^(/i, 2 ,f) would be solved by using equations (4.1), 
(4.2a) and (4.2e). Consider the Laplace equation (4.1). In this order 0(e 2 ), there are 
nonlinear forcing terms involved. In the process of solving $^(/i, 2 , t), it is recognized 
that an additional potential term without 2 -dependence, denoted by ^m(/*, 0> should 
be added to balance the system. Specifically, the form of the nonlinear forcing terms 
in the kinematic condition together with the requirements of the conservation of mass 
condition requires an additional contribution to the function $^(//, 2 , t). 
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Therefore we may assume that $ (2) (/z, 2 , <) to have the following form. 




It is remarked that <£$(alO (as well as <t> {2) ) must satisfy any conditions required 
on z,t). 

Thus (4.1) becomes 

VV (2) (^,z,0 = 0 (4.3) 

V 2 ^(/i,0 = 0 (4-4) 

Equation (4.3) can be solved, in a similar way to the case of 0(e), with the frequency 
of t-dependence being “2”. 


<j>W(fi,z,t) = sin(2t) 7 m/o(fm/*)cos/ m (z + A), (0 < /1 < 1, -A < z < A). (4.5) 

m=0 


To solve (4.4), we let (/*,<) = T(t)i jj’f/i). Then (4.4) yields 

Pi'S . Igjff „ 

dn 2 /i dp 


(4-6) 


Integrate both sides of (4.6) to obtain 


which implies that 



= — Infi + Const., 



— —Const. 
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By integrating both sides of the equality above, we obtain the solution of (4.6): 


4> { m = In/*, 


(4.7) 


where E\ is a constant. Set T(t) = sin2i. Then the solution of (4.4) is 


4>m ] = sin(2t)£’i In/*, 


(4.8) 


Physically, this is a source term correction to $^(/*,z, t). Therefore, the solution 
form for is: 

0 = sin(2<) | 7 m / 0 (/ m ^) cos l m (z + A) + f?i In /ij . (4.9) 

lm=0 J 

Considering the form of the nonlinear terms of 0(e 2 ), in the kinematic and normal 
force conditions, it may be assumed that the appropriate form of the solution for F ^ 
has a time-dependent term and a steady state (time-independent) term. The time- 
independent term will balance certain the nonlinear terms in the normal force balance 
at this order. Thus, by using similar techniques in solving this problem as were used 
in the linear order, the appropriate solution form for F ^ is assumed as follows: 


'1 F <21 (r,() = cos(21) <mCOs(„(z + A) + £ < m cos/ m (z + A). 


(4.10) 


m =0 


m=0 


Let 


F^ x \z,t) = cos tf^ x \z) and $^(/i,z, <) = sint$^(/*,z). (4-11) 
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Substitute these-expressions in (4.11) into the kinematic condition (4.2b) to get 


0$(2) m - m - m d 2 4>< l > • n d¥ l) dF^ 

2w (0) sin2<F (2) + sin2f-^- = -2u.’ (1) sin tF (l) -sm'2tF { ) + sm2t—_ 

(4.12) 

Substitute these expressions in (4.11) into the normal force equation (4.2c) and 
then partially differentiate both sides with respect to t to get 

d 2 F (2) 


— 4u^°* sin 2t¥ 2 ^ + sin 2 tF {2) + 2sin2f- 


dz 2 


d¥v . n fd¥ l) \ 


= 2u;^ sin + 2u/ 0 ^ sin 2tF (,) — r sin2f 

OH 


\W) 


— sin 2 1 


m 


+ sin 2f(F^) 2 — sin 2 1 


(f) : 


(4.13) 


Combine (4.12) and (4.13) to get 

a*( j ) d 2 (d¥ 2) \ 

-4w' 0 > sin 2i*W - sin 2t— - sin 21^ (-gjp j 


_,,#$(!) a4 (1) aF (1) 

= 2w^ sinf$^ + sin2fF^ — rr-z sin2F 


dfi 2 


dz dz 


+2n,<»sint^+sio2<(^") 3 - sin2(^j 


92 (d¥ l) dF^\ 

+am2t d? [- dT-dT j 

+2u> (1) sin 1 4 (1) + 2u> (0) sin 2*F (1) - sin 2 1 
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-sin 2 1 + sin2<(F (1) ) 2 -sin2< 


(4.14) 


Rewriting (4.14) gives 


0 = sin tu> 


(i) 


4$ (1) + 2 


d 2 F {l) 


dz 2 


. rt ( , 0) - (2) d$ (2) d 2 [d*W\ 

+ sin 2t |4u; $ + Qfl + q z i [ dp ) 


fi*) 


dF (1) 


d^ 2 


dz dz 








-(¥)! > 


(4.15) 


By (4.15), in order to vanish the secular term, the coeficient of sin t has to be zero. 
Since 

,83 fO) 

( 4 . 16 ) 


4 *<" + 2 ^^- * 0 , 


we must have 


dz 2 


J" = 0. 


(4.17) 


Applying the triple contact line condition (4.2g) for the order of 0(c 2 ), and by 
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(4.18) 


(4.10), orie obtains, when z = A 


F {2) { A, f ) = cos(2f) Y S m cos mir + ^ cos mi = 0, 


m=0 


m =0 


and when z = — A 


OO OO 


F (a) (-A,0 = COS(20 Y + E = 0. 


1=0 m=0 


(4.19) 


Combine (4.18) and (4.19) to get 


OO OO 


cos(2<) ^ 2fc — 

fc=0 t=o 


(4.20) 


Subtract (4.19) from (4.18) to get 


cos(2f) Y + Y = ®- 

k=i k= l 


(4.21) 


Since (4.20) and (4.21) are valid for all t, we must have the following systems. 


and 


Y ^ k 

t= o 


Y ^ 2k ~ i 

y k = 1 


Y 

k—0 


Y 

y k=i 


= 0 


= 0 , 


= 0 


= 0 . 


(4.22a) 


(4.226) 
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By the conservation of mass condition (4.2f), and by (4.10) 


f 

J - A 


cos(2t) 53 6 m COS l m (z + A) 4- 53 6 m cosl m (z + A) + C °^ — (F (1) ) 


m =0 


m=0 


+ 




dz = 0, 


which is valid for all t. It follows that both 


and 


/ [£ 6 m cos l m {z + A) + ]- 

■ / - A Lm=0 L 


/ A [f^ m cos/ m (z + A) + i(F (l >) : 

J — A n ^ 


dz = 0, 


m=0 


dz = 0. 


Hence, we have the following conclusion. 


6 0 = 6o = ~ 


Apply the kinematic condition to obtain 


2u > (0) 53 COS / m (z + A) + £ COS / m (z + A) + F t 


msO 


m=0 


: _».a^ (11 _ + ai"> a/-'" 


& 


where at^ = 0. 


dpi 1 


dz dz 


(4.23) 


(4.24) 


Using the orthogonal properties to eliminate the z-dependence, one obtains 

02$(i) d$(») dFW 


£ ' = -^-2l/3h : 


dpi 1 


dz dz 


dz, 


(4.25) 
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( 4 . 26 ) 


and for m~=':l,2,' • •, 


\Lli(i m hm + 2J 0) \6, 


r\ #2a(1) 

= ■// "-gs-'°*U*+W: + J_ 


A d$ {l) dF {l) 

~dz dT cos /m(; + A)d: - 


It is noted that the use of the orthogonal properties involves multiplication through 
by cos l q {z + A) and integration over the range (-A, A), with the appropriate (integer) 

q value. 

By the normal force balance condition, we have 




dz’ 


d\idt \ dn ) 

_(a|W) _ 2 [fOf + (<*!!!) + B ,v<>) (427) 

For the normal force condition without {-dependence, the orthogonal properties 
are used to obtain 

B N m = 4 1 r >«■ 

2A J-k dfi 2 \ dp J 


* H?)' 


and for m = 1,2, 


6 m = 


A(1 


— f‘ 

-idj- 


A 

)y- A 


w (o)/K.)^ + lf^V 

dn + 2 V dtx ) 


(4.28) 
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-H?) 


cos l m (z + A )dz. 


(4.29) 


For the normal force condition with ^-dependence, the orthogonal properties are 
used to obtain 


7o = 


4w<°> 


— r 

<°>A J- A 


0/* 2 v a? ) 




dz 


(4.30) 


and for m = 1,2, • 


2Aw*°*/o(/ m )7m - A(1 - £)«„ 


/ A 1 M /fliO >\ J 

“ L F cos ,m(2 + A)<b + 5 L (“arj cos '" (z + A),i!r 

i /a ( d &')\ 2 

+ 2 /- A (“«rj 00,< - (,+A)<£ * 

— / a(^ < 1) ) 2 c °s/ m (^ -h A)</z + ^ / A (^7~) COs/ m( Z + A)<fe. (4.31) 

By the Lanzcos Tau method and using (4.31), (4.26) and (4.21a), a set of linear 
algebraic nonhomogeneous equations in and 7 m are developed. Using techniques 
of numerical linear algebra, the truncated system can be solved for 7 m and S m for 
each m 6 {1,2, • • • , A/}. Using (4.29) and (4.21b), 6 m can be solved also. 

Knowing of the 7 m and 6 m can then be used to construct $* 2 *(/i, z, t) and F^(z, t). 
It is the shape function F ^ which is of most interest at second order, for it represents 
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the first nonlinear correction to the deformed interface shape of the finite length liquid 
column in natural harmonic oscillation. Also, the steady state correction to F 
indicates a modification of the mean form. 
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Chapter 5 Results at Second Order of 0(e 2 ): 
The Shape Function and Velocity Potential 

Corrections 


In this chapter, some numerical results at the second order of 0(e 2 ) are displayed. 
Results for the first six modes ( p = 2, 3, 4, 5, 6, 7) are presented. For each of these 
modes, the shape function F ^ in 0(e 2 ) is computed The initial parameters a* 1 * and 
0^ are so chosen that either = 1 and (3^ = 0, or = 0 and 0^ = 1. 
The value of e is set to be 0.4. With these values of the parameters, the corrected 
deformation up to 0(e 2 ) is plotted. 

Note that 

= F<°\z,t) + c FW(z,t) + jF«H*,t) 

= l+eFf'Hz,t)+‘^F»\z,t) (5.1) 

The perturbation contribution to F from both order 0(e) and 0(e 2 ) is graphed 
for various modes and A values. 
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FREE SURFACE PROFILES: LINEAR AND SECOND ORDER CORRECTED 

lambda=2.4, alpha1=0, beta 1 =1.0, p=4(mode), time=2*pi 
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FREE SURFACE PROFILES: LINEAR AND SECOND ORDER CORRECTED 

lambda=2.4 l alpha1=0, betal = 1.0, p=6(mode), time=2*pi 
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Chapter 6 Derivation of the Solvability 
Condition at Order of 0(e 3 ) 


The nonlinear correction to the interface shape has been determined at order t , 
and the forms of the theoretical solutions have been presented. 

However, it remains to investigate the nonlinear corrections to the interface for 
various families of parameters. 

Of interest is how the shape is modified by nonlinear corrections for various values 
of the slenderness parameter, for an even or odd number of half-wave deformations 
(at the linear order), for higher order modes in general, and for various forms of the 

initial disturbance (a (l) and /? (1) values). 

The third order (in e) system has been listed in Chapter 2. It is repeated below. 
In this order, the corrections to the time frequency, will be solved. 

The system equations are: 
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The Laplace equation: 


V 2 $ (3) (/i,z,t) = 0, (0 < n < 1, -A < z < A). 


Radial velocity condition: 


d$< 3 > 

dn 


= 0, {n = 0, 


-A < s < A). 


( 6 . 1 ) 


(6.2a) 


Kinematic condition: 


-J°> 


dF& 5$(3) 


dt 


+ 


dfi 




dt 

_3(F<'»)^ + 3 


dt ” 5/i 2 

d$ (2) d^O) 

3z dz 


<V 




dz d/idz 


dz dz 


(6.2b) 


Normal force equation: 


u> 


( 0 ) 


£$(3) 

dt 


— F( 3) _ 


d 2 F& 

d? 




dt 


dt 


dftdt 


-3^(0) F d)^^ _ 3u ,<°)/*3>. d2 * (1) 


dfidt 

-3u; (0) (F (1) ) 3 ??S -3 


dfidt 

d$W d$( 2 ) 
dpPdt '* dp dfi 


6jF m d* (l) d 2 * (>) _ 3 ^ (1) d$ (2) 


cfyi 5/i 2 


dz dz 


-6F (l) ^|^^^ -6F (1) F (2) 
az a/ioz 
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(6.2c) 


BF^V BF^BF™ 


+W _ 3f <»(^) + 3 


Bz Bz 


_ 9 (^)^ + e , m , , = 1 , 


Periodicity in time: 


V$ (3) (/r,z,t + 2 t) = V*(/i,z,t). 


Normal velocity: 


^$( 3 ) 


Bz 


= 0 . 


z=±\ 


Conservation of mass: 


J X |^F (3) + F (1) F (J) j dz = 0. 


Triple contact line condition: 


F (3) (±A,t) = 0. 


Assume the appropriate solution forms of $ (3) and F (3) as 


$ (3) (/i, Z,f) = 7m ) (0M^«nM)cOs/ TO (z + A), 

m=0 

F< 3 >(*,0 = £ 4 3) (0cos/ m (z + A). 


msO 


This is consisted with all boundary conditions. 


( 6-240 

(6.2e) 

( 6 . 2 /) 

(6.2«) 

(6.3) 

(6.4) 
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Let 


— 


<t> ( 1 ) (/i,z,t) = sin t 4> (1, (//, z), 


— 


$ ( 2 ) (/i,z,f) = sin 2 1 [4> (2) (/x, z) + Ei In /i] , 


— 


F ( 1 ) (z,t) = cos tF {l \z), 


— 


- ( 2 ) 

F w {z,t) = cos 2 tF {2) {z) + F ( 2 ). 


— 

Then substitute (6.3) and (6.4) into the Kinematic condition ( 6 . 2 b) to get 

— 


E ^r^co si m (z + a) + E ySHt)UK(L)cosi m (z 

m =0 aC m =0 

+ A) 

— 

= 

sintf/ftfj] + sin 3 i[ A' ^ 3 ] T 

(6.5) 


where 




[KH,\ 

= - ^1 + 

2 [ dn 2 l J ' 2 dn 2 




xwa 2 ** 1 ) 3 3 afw a$< 2 > 

5/i 2 4 * ' dp 3 + 2 dz dz 


— 


* ( 2 ) . 

ZdFWdQW n dF 
4* fW + 3 

2 5z 5/i3z 2 3 z 5z 5z 5z ’ 

( 6 . 6 ) 


and where 





[KH 3 ] 

“ 2 [ d/i 2 J 2 a^ 2 4 ^ J a/ 1 3 





3 a^ 1 ) a$< 2 > 3 dF^ a 2 #* 1 * 3 aF< 2 > a** 1 * 

2 dz dz 2 az a/iaz 2 az az 

(6.7) 
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Substitute (6.3) and (6.4) into the normal force balance equation (6.2c) to get 




m=0 


dt 


■/„(/») cos U* + A) - £(1 - O&0 co»U* + A) 


771 =0 


= cos tlNHi] + cos 3 t[I<H 3 ] + B\ {3) 


( 6 . 8 ) 


where 


= -3w (2) 4> (1) -3 u/ (0) F (,) 


(fH -I--'"? 

5$(2) 


3 - ( 1 ) d$W d 2 $W _ 3c^d$^ _ 3 -jn dfr^a 2 ^ 1 ) 

2 d/i d/i 2 2 dz dz 2 dz dpdz 

-3 F (l) F< 2) - 6 F ( 1 ) F 2) + ^(F (1) ) 3 - ^F (1) (^ 7 -) 

* (2) „ Q 

3 dF (1) dF (2) dF (1) dF 27 / dF (1) \ d 2 F (1) 

^2 dz dz dz dz 4 \ dz / dz 2 


(6.9) 


and where 

[JVtf 3 ] = 


-3u> (0) F (1) 


3 d4 (1) 
+ 2 d/i 


d4< 2 > 


[ dp 

d£( 2 ) 


+ £,] - 


d/i 


+ E\ 


3 d$W d 2 ** 1 * 3 d$< x > d$< 2 > 

2 d/i d/i 2 2 dz dz 


"^2 dz d/idz 


_ + 1 ( #<.))3 _ 3p ( i) p|yiy 


+ 


3 dF< x > dF< 2 > 9 / dF (1) \ 2 d 2 F<*> 

2 dz dz 4 \ dz / dz 2 


( 6 . 10 ) 
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— 

Taking derivative with respect to t in the normal force balance equation, we 

obtain 

— 

E d 2 " l '£ {t) Io(L)cost m (z + \)- £ (‘ - cos W.- + A) 

m-Q m = 0 


— 

dB.\ {3) 

— — sin tfiV/Zj] — 3 sin Zt[K H^\ 4- ^ 

(6.11) 

— 

Apply the orthogonal properties to the kinematic condition and to the normal 


force balance equation to eliminate the r-dependence. 



Thus the kinematic condition becames the following: for m = 0, 


— 

— 2A w (°>^ 3 i >( * ) + 0 = sint /' [KHi]dz + sin3f 1* [I<H 3 }dz, 

dt J- A J- A 

(6.12) 


and for m > 1 , 


— 



— 

= sin tf [KH\}cosl m (z + \)dz + sin3< / [K H 3 ] cos l m (z + \)dz. 
J-\ 

(6.13) 

— 

And the normal force balance equation becomes the following: for m = 0, 


— 

2A^«-2A 

at 1 dt 


— 

,K M , /A dBN^ , 

= -sin tj JNHxWz - 3 sin 3 1 J SH 3 ]dz + J ^ —^—dz. 

(6.14) 

— 

and for m > 1, 


— 



— 

= —sin t f [N H l ) cos l m {z + \)dz- 3 s\n3t f [NHs]cos l m (z + \)dz. 
J- A ■'“A 

(6.15) 

— 
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— 

0*01 NAL PAG* ft 

or poo* QUAunr 



Note that equation (6.15) can be rewritten as 

■#£»(<) , : „./ A ltftf.1 


dt 


— slot J 


»A(1 -PJ 


cos l m (z + \)d: 


, , . / A IJVHj] 

+3 sin 3t / — —cos /, 

y-A A(i - ft) 


, , . w . . (01 /0 (lm)<P^(t) 

,{z + d0 


Substituting normal force banlance equation into kinematic condition, we 
the following: for m = 0, 

J 3 )/ 


-2Au,<°> 


**•/>•>* 


3 /A 1 M 

+ 2A - 2A L ~ar iz 

= sint / A + sin3t / A [/C//^- 

«/ — A 


This can be rewritten as 


-2A(u.< 0 ') 1 


*iPW 


dt J 


= sint / A {[A:^] + - (0) (^i]}^ + S in3t {[KH 3 ] + 3u^[NH 3 ]} dz. 


For m > 1, 

-Au,<°> 




+3 sin 3t / - v cos l m (z + \)dz 

•/-A A(1 - ft) 


+ A/ ro /'(/ m bi 3) (0 


= sin 


tf [KHi] cosl m (z + A)dz + sin3t f [K Hz] cos l m (z + A)dz. 

«/ — A •/— A 


(6.16) 

obtain 


(6.17) 
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This can.be rewritten as 


—A + A/ m r o (/ m hS'(0 

<°>[A r //,] 


= sin tj^ [KHi\+ ^_ [2 


cos l m (z + A)Jz 


+ sin3t 


i: 


[KH 3 ] + 


3w (0) [A r // 3 ] 

1 -/ 2 , 


cos / m (z 4- A)dz. 


(6.18) 


Equation (6.18) is an ordinary differential equation for yj£K 
In order to get the solvability condition, set the coefficients of sin t to be zero. 
Therefore, for m = 0, 


f + J 0| (A'//,]} cos(0)<fc = 0, 


(6.19) 


and for m > 1, 


r\ ( u,(o) 'l 

/-a (1**1 + rrjrl^l} cos l m (z + A)dz = 0, m = 1,2,3, • • • (6.20) 


mr 


Note that / 0 = 0 and l m = Hence 

** i\ 


L{ 


o,(°) 1 

[KH\\ + cos l m (z + \)dz = 0, m- 1,2,3, 


1-^ 


( 6 . 21 ) 


By (6.6) and (6.9), the definitions of [KH\\ and [NHi], the solvability condition 
(6.21) can be rewritten as: 


0 = 


d 2 fr l) 

dpt , 2 
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Vf.ft£ 2 * (,) 3 ^ - (1) 2 5 3 4 > ^ 1) 3 dF<‘>d4>< 2 > 

_3F dy 1 4 ^ ’ dy 3 + 2 dz dz 

(2) . 

3 - m dF (1) d 2 4> (1) 3 dF< 2) d$ (1 > dF d$ (1) 

"*" 2 F dz d//dz 2 dz dz F dz dz 


+ 


u> 


(0) 


1 - 


-3u.«* l, > - 3u,l°>F l '> + E, 


) 


‘ 

2 dfi dfi 4 a/i^ 




a$(»a2^(i) 3 a^ (l) d^ <2 > 

~d~y d^ 7 " “ 2 dz dz 


d) ^ (1) d 2< E ( \ _ 3f’(i) _ 6F (1) F <2) + -(F (1) ) 3 
2 dz d/idz 2 V 


/’(i) 

4 


/d£W\ 3 3 dF (1) dF< 2) 

^ dz J **" 2 dz dz 


+3 


dF<*> dF 


( 2 ) 


dz dz 
It follows that 


/ dF (1)> 

2 d 2 F (1) 


V * , 

/ dz 2 

i 


> cos l m (z + A)<fz. 


( 6 . 22 ) 


( 2 ) = 


US' ’ = 




cos/ m (z + A )dz 
- 


£/-*.{[-! 


pW 


/d 2 $< 2) 

V“V~ 


d 2 ^ 1 ) ^(2)d 2 $< l > 3, ftlllx2 0 s * (I) . 3 dF (1) d4 (2) 


d/x 2 


d/i 3 


+ X 


2 dz dz 


-W - , , 

3 ^dF^d 2 ** 1 * 3dF (2) d<& (1) , 0 dF d$ (1) 

j__ FK 1 ) |- 3 — — 

2 dz dydz 2 dz dz dz dz 
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OMQINAL PAGC IF 
OF POOR QUALITY 



:+ 


jjj 


( 0 )- 


1 - ‘l 


-- W (0) F< 2) 

2 dfi 


-3 J 0) F {X) 

di> {l) 


'dfrv 

r^~ 


+f,) 


_ 3^<° >F - -u,<°>(F (, >) 2 V— 

dji 4 a/i 2 

3 di>M (d& 2) \ 3 - n) d¥" 3 d$< 2 > 

~2 dp { dfi JrEl )~2 F dn dn 2 2 dz dz 

--F (1) -| - (l) ^ - ^ - 3F (1) F (2) - 6F (l) F 2) + ^(F (,) ) 3 

2 az OfiOz 2 


-^(i) 

4 


+3- 


(0 


3 dF (1) dF (2) 
2 #2 dz 


x(2) 


/dF< l > N 

l 2 Qlp( 1) 


V dz j 

I dz 2 



$2 C?2 4 

Therefore, cj^ can be determined numerically. 


> cos l m (z + A )dz. 
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OftiaiNAL PA3€ fi 
OF rOO« QUALITY 


(6.23) 



Chapter 7 Results of Third Order of e, 0(e 3 ): 
Corrections to the Frequency, u 


Preliminary calculations have yielded the corrections to u i. Note that 

UJ = jo) + t l u m. (7.1) 

Results are plotted for — (0) — versus the slenderness parameter A. This is done 
for modes p = 2,3 and 6. 


52 



RATIO OF NONLINEAR FREQUENCY CORRECTION TO LINEAR FREQUENCY VS LAMBDA 
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RATIO OF NONLINEAR FREQUENCY CORRECTION TO LINEAR FREQUENCY VS LAMBDA 

P-6 (mod*), total -1.0 
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Appendix A 

Code for Linear Root Finding 
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* 


* 


* 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* 


* * 


* * * 


* * 


FROM THE FIRST ORDER OF EPSILON WE HAVE TWO EQUATIONS WHICH 
GIVE TH&- RELATIONSHIP BETWEEN THE SLENDERNESS LAMBDA AND THE 
ANGULAR FREQUENCY OMEGA (0) AT THE ZERO ORDER OF EPSILON. 

ONE OF THE 'EQUATION IS FOR EVEN MODE. THE OTHER ONE IS FOR 
ODD MODE . 

HERE WE CALCULATE THE ODD MODE EQAUTION. 
************ 


DIMENSION BL(80) , B1 (80) , BO (80) 

REAL* 8 BL , B1 , BO , A , PI , LAMBDA , S , S 1 , S AVG , 
* X , XMAX , DELTX , FMAX , XAVG , FAVG , XI , FX1 , FX 
INTEGER N , I , K 


C 

C INITIAL THE VALUE OF LAMBDA 

C 

LAMBDA=1. 80 DO 
WRITE (6,1) LAMBDA 
1 FORMAT ( ' ' , ' LAMBDA= ' , F7 . 4 ) 

PI=3. 1415926D0 

A=(0.5D0) *( (LAMBDA-DTAN( LAMBDA) )/DTAN( LAMBDA) ) 

C 

£************ 

C CALCULATE THE WAVELENGTH BL(N) , THE MODIFIED BESSEL FUNCTIONS 

C OF ZERO AND FIRST ORDER OF FIRST KIND I0(X) AND I1(X). 

c ************ 

C 

DO 20 K=l, 80 

BL(K)=(K*PI)/(2.D0*LAMBDA) 

BO(K)=BESSIO(BL(K) ) 

Bl(K) ss BESSIl (BL(K) ) 

20 CONTINUE 



200 


* * * * 


******* 


USE BISECTION METHOD TO SOLVE THE ODD MODES EQUATION. 

DEFINE FX TO BE THE LEFT-HAND SIDE OF THE EQUATION. 

OMEGA IS DENOTIED BY X WHICH WILL BE THE ROOTS OF THE EQUATION. 


THE INITIAL INTERVAL WILL BE CHOSEN AS [X, XMAX] . 

DELTX IS CHOSEN AS 0.005. 

LET N=80 BE THE NUMBER OF BISECTIONS DESIRED. 

CONSIDERING THE FX GOES TO INFINIT WHEN FX > FMAX (-10000) . 

***** ****** 

READ(5,200) X,XMAX 

DELTX-0.005 

FMAX=10000 

N=80 

FORMAT (4F10.4,I2) 


C 

C 

C 

C 

C 

C 

C 

C 


************ 

FROM THE LEFT END VALUE OF THE INITIAL INTERVAL, X, WE ARE 

GOING TO HAVE THE FUNCTION FX. 

THE LOOP IS TO ADD 40 TERMS FOR THE SUMATION WHICH IS 
DEFINED AS S. 

************ 
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0000000000(1000000 


S=0 . DO 

DO 201 K= 1,20 

S=S+ ( 1 . DO/ Cl • D0-BL(2*K) **2 ) ) * ( 1 • DO/ (X**2+BL(2*K) * 

* (l.D0-BL(2^K) **2) *(B1(2*K)/B0(2*K) ) ) ) 

201 CONTINUE 

300 FX=A-S* (X**2) 

C 

c ************ 

C NOW WE DEFINE A NEW VALUE AS X1=X+DELTX. 

C REPLACE X BY XI FOR FX. THEN WE HAVE FX1. 

C THE LOOP IS TO ADD 40 TERMS OF THE SUMATION, WHICH IS 

C REDEFINED AS SI. 

C************ 

C 

400 X1=X+DELTX 
S1=0. DO 

DO 401 K= 1,20 

S1=S1+ ( 1 . DO/ ( 1 . D0-BL(2*K) **2) ) * ( 1 . DO/ (X1**2+BL(2*K) * 

* ( 1 . D0-BL(2*K) **2)*(B1(2*K)/B0(2*K) ) )) 

401 CONTINUE 
FX1=A-S1*(X1**2) 

********* *** 

NOW WE NEED TO CONSIDER THE TWO RESULTS FX AND FX1. 

IF FX*FX1 IS LESS THAN ZERO, THE ROOT MUST BE IN THE INTERVAL 
[X, XI]. (WHICH MEANS FX AND FX1 HAVE DIFFERENT SIGNS.) 

IF FX*FX1*0, THEN XI IS THE ROOT OF THE FUNCTION. 

IF FX*FX1 IS GREATER THAN AERO, THE ROOT MUST BE IN THE 
INTERVAL [XI, XMAX] . 

FURTHER, IF FX*FX1<0, LET XAVG= (X+Xl)/2 . 

IF FX*FX1>0, LET XAVG= (Xl+XMAX) /2 . 

BY SUBSTITUTING XAVG WE WILL GET FAVG. (DEFINE AS SAVG 
FOR THE SUM OF THE FIRST 40 TERMS OF THE SUMATION.) 
CONTINUE THE SAME PROCEDURE UNTIL THE 'REAL ROOT',X, 

IS OBTAINED WHICH MAKES THE FUNCTION FX ZERO. 
************ 

IF (FX*FX1) 800,500,700 
500 WRITE (6, 600) XI 

600 FORMAT ( ' ' , • X* • , F24 . 18 , 1 IS A REAL ROOT') 

X-X1+DELTX 
GO TO 300" 

700 IF (XI. GE. XMAX) STOP 
X=X1 
FX=FX1 
GO TO 400 

800 DO 1100 I-1,N 
XAVG= ( X+Xl ) / ( 2 . DO ) 

SAVG=0 . DO 
DO 801 K“1 , 20 

SAVG=SAVG+ ( 1 . DO/ ( 1 . DO-BL ( 2 *K) * *2 ) ) * ( 1 . DO/ (XAVG**2+ 

* BL( 2*K) * (l.D0-BL(2*K) **2) * (B1 (2*K) /BO (2*K) ) ) ) 

801 CONTINUE 
FAVG=A-SAVG* (XAVG**2) 

IF(ABS(FAVG) .GT.FMAX) GO TO 1400 
IF(FX*FAVG) 1000,1200,900 

900 X=XAVG 
FX=FAVG 
GO TO 1100 
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1000 X1=XAVG 
FX1=FAVG 

1100 CONTINUE- - 

1200 WRITE (6 jS00X~‘XAVG 
1300 FX=FX1 - 
X=X1 

GO TO 400 

1400 WRITE(6, 1500) XAVG 
1500 FORMAT ( ' 'FUNCTION 
GO TO 1300 
END 


APPROACHING 


INFINITY FOR X= 


, F7 . 


C 


FUNCTION BESSIO (X) 

RETRURNS THE MODIFIED BESSEL 


10 FOR ANY REAL X. 


C 


C 


C 

C 


C 

C 

C 


C 


REAL* 8 Y, PI, P2, P3 , P4 , P5, P6 , P7 

accumulate polynomials in doulble PRECISION 
REAL* 8 Ql, Q2 , Q3, Q4, Q5, Q6, Q7 , Q8 , Q9 


P1=1.0D0 
P2=3 . 5156229D0 
P3=3 . 0899424D0 
P4=l . 2 067 49 2 DO 
P5— 0 .2659732 DO 
P6=0. 360768D-1 
P7-0.45813D-2 


Q1=0 . 39894228D0 
Q2-0.1328592D-1 
Q3=0. 225319D-2 
Q4—0.157565D-2 
Q5-0.916281D-2 
Q6— 0.2057706D-1 
Q7=0 . 2635537D-1 
Q8— 0.1647633D-1 
Q9=0. 392377D-2 


THEN 


POLYNOMIAL FIT 
[F ( DABS (X) . LT . 3 . 7 5D0 ) 

(P3+V* (P4 + Y*(P5«* (P6 + Y*P7) ) ) )) 

ELSE 

kX-DABS(X) 

BESSIO-%MP(AX)/DSQRT(AX) ) * (Ql+Y* (Q2+Y* (Q3+Y* (Q4 

+Y* (Q5+Y* (Q6+Y* (Q7+Y* (Q8+Y*Q9 )))))))) 


ENDIF 

RETURN 

END 


StURNS N THE S MODIFIED BESSEL II FOR ANY REAL X. 


REAL* 8 Y, PI, P2, P3 , P4 , P5, P6, P7 
ACCUMUIATE X POLYNOMIALS IN DOUBLE PRECISION 

REAL^8 Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9 

P1=0.5D0 


P2»0.87890594D0 
P3*0 . 51498869D0 
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P4=0. 15084934D0 
P5=0 . 2 6587 3 3D- 1 
P6=0.301-5-32B»2 
P7=0 . 3 2 4*1 ID- 3**- 

Q1=0 .39894228 DO 
Q2=-0 . 3988024 D-l 
Q3=-0 . 362018D-2 
Q4=0. 163801D-2 
Q5=-0. 1031555D-1 
Q6=0 . 2282967D-1 
Q7=-0 . 28953 12 D-l 
Q8=0 . 1787 654 D—l 
Q9=-0 . 420059D-2 

POLYNOMIAL FIT 

IF ( DABS ( X ) . LT . 3 . 7 5 DO ) THEN 

Y= (X/3 . 7 5 DO) **2 

BESSI1=P1+Y* (P2+Y* (P3+Y* (P4+Y* (P5+Y* (P6+Y*P7) ) ) ) ) 

BESS I 1=X* BESS II 

ELSE 

AX=DABS ( X ) 

Y=3 . 75D0/AX 

BESSI1=*(DEXP(AX) /DSQRT(AX) ) * (Ql+Y* (Q2+Y* (Q3+Y* (Q4 
*+Y* (Q5+Y* (Q6+Y* (Q7+Y* (Q8+Y*Q9) ))))))) 

ENDIF 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 


FROM THE FIRST ORDER OF EPSILON WE HAVE TWO EQUATIONS WHICH 
GIVE THE -RELATIONSHIP BETWEEN THE SLENDERNESS LAMBDA AND THE 
ANGULAR FREQUENCY OMEGA (0) AT THE ZERO ORDER OF EPSILON. 

ONE OF THE EQUATION IS FOR EVEN MODE. ANOTHER ONE IS FOR ODD 

MODE. 

HERE WE CALCULATE THE EVEN MODE EQAUTION. 


DIMENSION BL(80) ,B1(80) ,B0(80) 

REAL* 8 BL , B1 , BO , A , PI , LAMBDA , S , S 1 , SAVG , 
* X , XMAX , DELTX , FMAX , XAVG , FAVG , XI , FX1 , FX 
INTEGER N,I,K 


*** INITIAL THE VALUE OF LAMBDA 


LAMBDA=1.80D0 

WRITE (6,1) LAMBDA 

FORMAT ( ' 1 , ' LAMBDA* ' , F7 . 4 ) 

PI=3 . 1415926D0 

A=(0. 5D0)* (LAMBDA) *(DTAN( LAMBDA) ) 
************ 


* 


CALCULATE THE WAVELENGTH BL(N) , 
OF ZERO AND FIRTS ORDER OF FIRST 
******** 


THE MODIFIED BESSEL FUNCTIONS 
KIND I0(X) AND I1(X) . 

* * * 


DO 20 K=l, 80 

BL(K) * (K*PI) / ( 2 . D0*LAMBDA) 
BO (K) =BESSI0 (BL(K) ) 
B1(K)=BESSI1(BL(K) ) 

20 CONTINUE 


READ (5, 200) X, XMAX 

DELTX=0.005 

FMAX=10000 


200 


201 

300 

400 


401 

500 

600 

700 


N=80 

FORMAT (4F10 .4,12) 


S=0.D0 

s2s+(1.DQZ(1?DO-BL(2*K-1)**2) )*(1.D0/(X**2+BL(2*K-1)* 

* ( 1 . D0-BL(2*K-1) **2) *(B1(2*K-1)/B0(2*K-1) ) ) ) 

CONTINUE 
FX=A+S * ( X* * 2 ) 

X1=X+DELTX 


S 1=0. DO 

DO 401 K=l, 20 

S1=S1+ ( 1 . DO/ ( 1 . D0-BL(2*K- 

* (l.DO-BL(2*K-l)**2)*(Bl( 


1) **2) ) *(1.D0/(X1**2+BL(2*K-1) * 
2*K-1)/B0(2*K-1) ) ) ) 


CONTINUE 


FX1=A+S1*(X1**2) 

IF (FX*FX1) 800,500,700 
WRITE (6, 600) XI 

FORMAT ( 1 ' , 'X=‘ ,F24.18, • IS A REAL ROOT*) 

X-X1+DELTX 

GO TO 300 

IF (XI. GE. XMAX) STOP 


X=X1 
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800 


FX=FX1 
GO TO 400 
DO 1100 X=1 rJJ-. 

XAVG= ( X+Xl ) / C2 . DO) 

SAVG=0 . DO : 

DO 801 K= 1,20 

SAVG=SAVG+ ( 1 . DO/ ( 1 . D0-BL( 2*K-1 ) **2 ) ) * ( 1 . DO/ (XAVG**2+ 

* BL(2*K-1) * (l.D0-BL(2*K-l) **2) *(B1(2*K-1)/B0(2*K-1) ) ) ) 
801 CONTINUE 

FAVG=A+SAVG* (XAVG**2) 

I F ( ABS ( FAVG ) . GT . FMAX ) GO TO 1400 
IF ( FX*FAVG) 1000,1200,900 
900 X=XAVG 
FX=FAVG 
GO TO 1100 
1000 X1=XAVG 
FX1=FAVG 
1100 CONTINUE 
1200 WRITE (6,600) XAVG 
1300 FX=FX1 
X=X1 

GO TO 400 

1400 WRITE(6, 1500) XAVG 

1500 FORMAT ( ’ FUNCTION APPROACHING INFINITY FOR X=',F7.4) 

GO TO 1300 
END 

FUNCTION BESSIO (X) 

RETRURNS THE MODIFIED BESSEL 10 FOR ANY REAL X. 

REAL* 8 Y, PI, P2, P3 , P4 , P5, P6, P7 
REAL* 8 AX, X, BESSIO 

ACCUMULATE POLYNOMIALS IN DOULBLE PRECISION 
REAL* 8 Ql, Q2, Q3, Q4 , Q5, Q6, Q7 , Q8, Q9 
PI— 1. 0D0 
P2=3 . 5156229D0 
P3=3 . 0899424D0 
P4=l • 206749 2D0 
P5=0 . 26597 3 2 DO 
P6=0. 360768D-1 
P7=0.45813D-2 
C 

Q1=0. 39893228D0 
Q2=0.1328592D-1 
Q3«0.225319D-2 
Q4=—0 . 157565D-2 
Q5=0 . 916281D-2 
Q6*-0.2057706D-1 
Q7=0 . 2635537D-1 
Q8=-0 . 1647633D-1 
Q9=0 . 392377 D-2 
C 

C POLYNOMIAL FIT 

IF (DABS(X) .LT.3.75D0) THEN 
Y=(X/3 .75D0) **2 

BESSI0=P1+Y* (P2+Y* (P3+Y* (P4+Y* (P5+Y* (P6+Y*P7) ) ) ) ) 

ELSE 

AX=DABS (X) 
y=3 . 75D0/AX 

BESSIO* (DEXP (AX) /DSQRT (AX) ) * (Ql+Y* (Q2+Y* (Q3+Y* (Q4 
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*+Y*(Q5+Y* (Q6+Y* (Q7+Y* (Q8+Y*Q9) ))))))) 

ENDIF 

RETURN - -.V 

END 

FUNCTION BESSI1 (X) 

RETURNS THE MODIFIED BESSEL II FOR ANY REAL X. 

REAL* 8 Y, PI, P2 , P3 , P4 , P5 , P6, P7 
REAL* 8 AX , X , BESSI 1 

ACCUMULATE POLYNOMIALS IN DOUBLE PRECISION 

REAL* 8 Q1 , Q2 , Q3 , Q4 , Q5 , Q6 , Q7 , Q8 , Q9 

P1=0 . 5D0 

P2=0 . 87 890594 DO 

P3=0 . 51498869D0 

P4=0 . 15084934D0 

P5=0 . 2658733D-1 

P6=0 . 301532 D- 2 

P7=0. 32411D-3 

Q1=0 . 39894 2 2 8D0 
Q2=-0 . 3988024D-1 
Q3=“0. 362018D-2 
Q4=0. 163801D-2 
Q5=-0. 1031555D-1 
Q6-0.2282967D-1 
Q7=-0. 2895312D-1 
Q8-0.1787654D-1 
Q9— 0.420059D-2 

POLYNOMIAL FIT 

IF (DABS(X) .LT.3.75D0) THEN 

Y= (X/3 . 75D0) **2 

BESSI1=P1+Y* (P2+Y* (P3+Y* (P4+Y* (P5+Y* (P6+Y*P7) ) ) ) ) 

BESSI1=X*BESSI1 

ELSE 

AX=DABS (X) 

Y=3 . 75D0/AX 

BESSI 1= (DEXP (AX) /DSQRT (AX) ) * (Ql+Y* (Q2+Y* (Q3+Y* (Q4 
*+Y* (Q5+Y* (Q6+Y* (Q7+Y* (Q8+Y*Q9) ))))))) 

ENDIF 

RETURN 

END 
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Appendix B 


Solvability Condition: 
Alternative Formulation 
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Determining the Solvability Condition Formally 


Start with the basic governing equation at 0(e 3 ) which is 


V 2 $ (3) =0, on 0 < n < 1, 0 < 6 < 2ir, and — A < 2 < A. 
Multiply by $U) to integrate over the volume 

f = 0. 

J volume 

Use cylindrical coordinates to write dV = d[i(fid9)dz, and integrate in n 

^ 2 $( 3 ) 1 1 QPqW Q 2 <fr( 3 ) 


///..«{: 


+ 


+ 


Since 


dn 2 n dn n 2 dO 2 dz 2 

1 d 2 $< 3 > 

= 0, 


) 


dfi(fidO)dz = 0 


y 2 dO 2 

the basic governing equation can be rewritten as follows. 

:: j j } + II h m ^r ilxiMz 

^$(3) 


///„*<•»! 


dz 2 


dfidddz = 0. 


Use integration by parts to get the adjoint system. 
Denote 


= and y = 


dn 
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Then 


du = and dV = ^ 2$(3) 


d/i 


du 2 ' 


It follows from integration by parts that 

/•A r2ir rl lu d 2 $^ 

LI l 

ll 


= r / 2 "L<& (i) $h dedz- [ a r +p*w)dpd0dz 

J — A J 0 q J — A Jo Jo 


A 


= J J 




d0dz - / A I** f 1 $< 3) $ (1) d/id0dz 
„ d-A do do 


- / / / /i^^d/idfldz. 

J-A Jo Jo » » 

= LrL"*? 


dedz 


Jm=i 

ii 


- / a a £* | $ (i) * (3) d^dz + y A £* £ ^^d/idfldz 

- y A £* \ *< 1 >* (3) d0dz + y A £* £(^ ] + H^Q^dfidOdz 


Jm=i 

For the second integral, we have 

r A y 2 * y 1 w,^ 3) 

du 


/■a /-a* yi d$( 3 ) 

/ / / $ (1) ^r- d^id^di 

d-A do do £711 


= f p 




Tl 




rA r2ir rl d^ 1 ) 

-III ^ 3 — dfidOdz 

J- a Jo Jo a/i 


For the third integral, we have 

rA /*2ir /i > m <9 2 <$( 3 ) 


/ a r 2 * rl . % ^<T)V o ; 

/ / 1 ~T 2 dfidOdz 

•A Jo JO UZ l 

r2ir rl rA , v ^2*( 3 ) 

= / / / /i$ (1) dfidOdz 

do do d-A OZ* 
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= £'£ 


= f/ 

u<->Tl 

az 


1 A 

-A 

A 

J —A 


dfidd — <f>^ dzdfidO 


d^e 


-[/cVoK’* 13 ’ 

Let d?r = dft(fidd)dz . Then by 


A 

-A 


ray vgw 


5 2 $( X ) 




and by setting 


dz 

^$(3) 


= 0, 


z=±A 


dz 


= 0 , 


one obtains 

0 = /// 

J J J volume 

r A <94> (1) 


z=dhA 


/ A . rA f 2ir 


/A f2ir r 1 

+ f f * (3) 

J-A JO JO 


l ~diS~ 


■2«J, 1 >J dfidddz + J A ^j**[& l) *W) l o d0dz 

- c j: £ * {z) * { ?w d > + r £ 4 


dfidO 


.. -a: 


^$(3) 


5$(1) 

52 


r2ir r l rA 

dfidd+ / / / (i—~—dzdud$ 

Jo Jo J - a az 1 


d ' 2 $W 


J-A 

A / 2 *Lmd$ (3) 5$( x > 


5/x d\i 


$ (3) d0dz 


1,1=1 


/ A f2* /*1 v 

/ / $ (3) 
-A JO JO 

= Iff * (3) 

J J J volume 




5/i 2 


dz 2 


dfidOdz 


L,/$ (1) 
+ -* . 1} 


d\i 


2 ' ft M dz 2 


dV 
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dpL dp 


/ / *« 

V-A JO 


$( 3 ) 


<L6dz 


Jm=i 


By Laplace equation on 0(c) and on 0(e 3 ), 

V 2 $ (1) = 0 and V 2 $< 3 > = 0. 


Thus the solvability on p = 1 becomes 

r A t 2 * 


r r *<■»— - $' 3 »— 

J- a 7o dp 9fi \ 


dOdz = 0 


or 


27T 


r U,££!!l _ *<3,^1 

/-a I dll d/i 


dz = 0. 


Solvability on free surface is then 

jC 

Kinematic 


rA 

-A 


9/j dii 


dz = 0. 


ir,\dF^ d^ 3) 

-J 0) -^r~ + = sint [tftfi] + sin3t 


dt dp 


Normal force 


Ffl F&) 

w(0)! i?r _ F(3) ■ ~d^~ = cost i NH d + cos3 * [^^3] + bn (3) . 


* (3) = S COS (/ m (z + A)) + G(/i, f ). 

m=0 


Solvability condition 


f (3) = £ ^(t)cos(/ m (z + A)). 

m=0 


/ A a [$ (1) $<, 3 > - $^$ (3) ] dz = 0. 
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$ (1) = Y sin tA m (t)I 0 (l m fi) cos(/ m (z + A)). 


m=0 


$ (3) _ £ ^(^IoiU^COsiLiz + A)). 


m=0 


F (3) = Y, tSHfyoilml*) COS (l m {z + A)). 


m=0 


(C) 


(A) 


Y A m I 0 (lmf*)\v=l I COS (l m (z + \))^dz 

n J — A 

rn=0 

Y l mA m r 0 (lmH)\ li =l [ COs{l m (z + \))$ {3) dz = 0 

* J — A. 


yi / COs(/ m (z + A)) 

m=0 • y_A 

Y t m A m r 0 (l m ) / A cos(/ m (z + A ))* (3) dz = 0 

•'-A 


— u>^— — |-sin< [KHi] + sin3t [KH 3 ] 

at 


Y AmIo{lm) f cos(/ m (z + A)) [sin t [KHi] + sin3f [KH 3 ]} dz 

J —A 


m= 0 


00 M 

- £ J 0 'A m I 0 (L) J cos (/ m (z + A))<fz 


m=0 


- f: IM&L) r coa(U(z + A))*' 3 '^ = 0 

A */ — A 

m=0 

Apply normal force equation to get 

JO)?* p(3) _ ?_Y— = cos t[N Hi] + cos 3t[NH 3 ]. 

dt dz 2 
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